## Replication of the main analyses from:
##
## Christopher Adolph, Kenya Amano, Bree Bang-Jensen, 
## Nancy Fullman, Beatrice Magistro, Grace Reinke, Rachel Castellano,## Megan Erickson, and John Wilkerson, "The Pandemic Policy U-Turn:## The role of partisanship, public health, and race in decisions
## to ease COVID-19 social distancing policies in the U.S.", 
## forthcoming in Perspectives on Politics.
##
## This file:
## make Figures 2, 3, and 6 and computes the contents of Tables B1
##
## Estimate pooled stratified Cox proportional hazards of 
## the first indoor easing of state-level social distancing measures, 
## 16 April to 6 July 2020; and produce plots of hazard rates 
## and expected acceleration  
##
## Christopher Adolph
## 18 June 2021     faculty.washington.edu/cadolph

## Clear memory
rm(list=ls())

## Set random number seed for consistency in replication
set.seed(123456)

## Use fancy fonts?
fancy <- FALSE

## Load libraries
library(survival)      # For Cox ph and simulation
library(coxed)         # for AME predicted survival times
## note that a helper file is loaded to fix a typo in
## coxed2; see below
library(MASS)          # for mvrnorm()
library(formula.tools) # for lhs(), rhs()
library(tile)          # graphics; available from
## http://faculty.washington.edu/cadolph/software
library(RColorBrewer)  # Color selection
library(simcf)         # for extractdata; available from
## http://faculty.washington.edu/cadolph/software

## Load helper functions, including a fix for an apparent bug in coxed 0.3.3
##  (see helperAnalysis.R for more details)
source("helperAnalysis.R")

## Set up colors
col <- brewer.pal(7, "Set1")
nicered <- col[1]
niceblue <- col[2]
nicegreen <- col[3]
nicepurple <- col[4]
niceorange <- col[5]
nicebrown <- col[7]

## Load data
load("AnalysisData.Rdata", verbose=TRUE)

## Save new copy of data for additional processing for the baseline model
allsurvData <- mData

## Set earliest date = April 16 >=50 (check)
startDate <- 50
allsurvData <- allsurvData[(allsurvData$DaysSinceFCT>=startDate), ]

## Set end date = July 6 <=131 (check)
endDate <- 131
allsurvData <- allsurvData[(allsurvData$DaysSinceFCT<=endDate), ]

## Set source of COVID case, death, and test data
allsurvData <- allsurvData[allsurvData$COVIDSource=="JHU",]  # was JHU

## Select policy types
allsurvData <- allsurvData[(allsurvData$PolicyType=="GathRestrictRecomAny")|
                           (allsurvData$PolicyType=="AnyBusinessClose")|
                           (allsurvData$PolicyType=="RestaurantRestrict")|
                           (allsurvData$PolicyType=="BarRestrict")|
                           (allsurvData$PolicyType=="StayAtHome")
                          ,]

## Set up model specification, with clustering and stratification
model1a <- sdEasing ~
    DemGovernor +
    TrumpVote2016 +
    PctBlack +
    TestPositivityRateLagSlopeLM +
    logNewDeathsPerM7DaysMAVG +
    NoDeaths7DaysMAVG +
    NewConfirmedCasesPerMLagSlopeLM +
    log(PopDensity) +
    strata(PolicyType, SubstatePolicyDateSubNonReligIndoorEasing) +
    cluster(as.factor(StateName))

## Create a dataset with only at risk cases
obsSurvData <- allsurvData[!is.na(allsurvData$PolicySubNonReligIndoorEaseDV),]

## Create survival object
sdEasing <-  Surv(time=obsSurvData$DaysSinceFCT-1,
                  time2=obsSurvData$DaysSinceFCT,
                  event=obsSurvData$PolicySubNonReligIndoorEaseDV
)

## Renumber rows (required by coxed below)
rownames(obsSurvData) <- 1:nrow(obsSurvData)

## Estimate CPH model: this is the baseline model
res1a <- coxph(model1a, data=obsSurvData, x=TRUE, y=TRUE)

## This summary is used to fill in the goodness of fit for Table B1.
## However, the counterfactuals for covariates come from the
## simulations below for Figure 2, as they mostly involve continuous
## covariates, for which the one-unit-change hazard ratios are unhelpful
summary(res1a)
AIC(res1a)

###################################################################################
## Set up counterfactuals for Figure 2

## Coefficient numbers
xNum <- c(1, 2, 3, 4, 5, 7, 8)

## Scenario names
scenname <- c("Republican Governor",
              "More Trump Voters",
              "Higher Black Population",
              "Downward Trend Positivity",
              "Lower Rate of New Deaths",
              "Downward Trend New Cases",
              "Lower Population Density"
              )

## Note: the next two commands must deal with an issue in computing quantiles 
## (in particular, the 25th percentile) created by replacing 0s with 1s 
## in the unlogged death rate data:  We use the uncorrected origNewDeathsPerM7DaysMAVG
## to compute the correct quantiles, while preserve the replaced values in 
## logNewDeathsPerM7DaysMAVG for in-sample simulation via coxed
            
## "Before" scenario (theoretically discourages easing)
## Unlogged versions of the below are also used 
## to fill in the "pre" values in Table B1
xpre1 <- with(obsSurvData,
               c(1,
                 quantile(TrumpVote2016, probs=0.25),
                 quantile(PctBlack, probs=0.25),
                 quantile(TestPositivityRateLagSlopeLM, probs=0.75),
                 log(quantile(origNewDeathsPerM7DaysMAVG, probs=0.75)),
                 quantile(NewConfirmedCasesPerMLagSlopeLM, probs=0.75),
                 log(quantile(PopDensity, probs=0.75))
                 ))                

## "After" scenario (theoretically encourages easing)
## Unlogged versions of the below are also used 
## to fill in the "post" values in Table B1
xpost1 <- with(obsSurvData,
              c(0,
                quantile(TrumpVote2016, probs=0.75),
                quantile(PctBlack, probs=0.75),
                quantile(TestPositivityRateLagSlopeLM, probs=0.25),
                log(quantile(origNewDeathsPerM7DaysMAVG, probs=0.25)),
                quantile(NewConfirmedCasesPerMLagSlopeLM, probs=0.25),
                log(quantile(PopDensity, probs=0.25))
                ))

## Make Figure 2; which must then be polished in Adobe Illustrator.
## The object hr is also used to fill in Table B1 hazard rates 
## and confidence intervals, and to label points in the polished plots
hr <- plotCovidCPH("Figure2",  res1=res1a,
                   scenname=scenname, xNum=xNum, xpre1=xpre1, xpost1=xpost1, 
                   limits=c(0.1, 3.1),
                   at=c(0.5, 1, 1.5, 2, 2.5, 3.0),
                   labsX=c("0.50", "1.00","1.50", "2.00", "2.50", "3.00"),
                   labsT=list(c("half", "as", "likely"),
                              c("just", "as", "likely"),
                              c("50%", "more", "likely"),
                              c("100%", "more", "likely"),
                              c("150%", "more", "likely"),
                              c("200%", "more", "likely")),
                   sortBy=1, rev=TRUE, logaxis=FALSE,
                   titleT="The first move to ease social distancing is...",
                   titleX="Hazard Ratio",
                   pch=c(16, 16, 15, 17, 15, 16, 16),
                   col=c("black", niceblue, nicered, nicepurple, nicered, niceblue, niceblue),
                   fancy=fancy
                   )

###################################################################################
## Make Figure 3: Cox duration marginal effects

## "Before" scenario for coxed (conditions that discourage easing)
cf2 <- c("TrumpVote2016=quantile((allsurvData$TrumpVote2016), 0.25), DemGovernor = 1",
         "logNewDeathsPerM7DaysMAVG=log(quantile(obsSurvData$origNewDeathsPerM7DaysMAVG, probs=0.75)), NoDeaths7DaysMAVG=0, NewConfirmedCasesPerMLagSlopeLM=quantile(obsSurvData$NewConfirmedCasesPerMLagSlopeLM, probs=0.75), TestPositivityRateLagSlopeLM=quantile(obsSurvData$TestPositivityRateLagSlopeLM, probs=0.75)",
         "DemGovernor = 1",
         "TrumpVote2016=quantile((obsSurvData$TrumpVote2016), 0.25)",
         "PctBlack=quantile((obsSurvData$PctBlack), 0.25)",
         "NewConfirmedCasesPerMLagSlopeLM=quantile(obsSurvData$NewConfirmedCasesPerMLagSlopeLM, probs=0.75)",
         "logNewDeathsPerM7DaysMAVG=log(quantile(obsSurvData$origNewDeathsPerM7DaysMAVG, probs=0.75)), NoDeaths7DaysMAVG=0",
         "TestPositivityRateLagSlopeLM=quantile(obsSurvData$TestPositivityRateLagSlopeLM, probs=0.75)",
         "PopDensity=quantile(obsSurvData$PopDensity, probs=0.75)"
    )


## "After" scenario for coxed (conditions that encourage easing)
cf1 <- c("TrumpVote2016=quantile((allsurvData$TrumpVote2016), 0.75), DemGovernor = 0",
         "logNewDeathsPerM7DaysMAVG=log(quantile(obsSurvData$origNewDeathsPerM7DaysMAVG, probs=0.25)), NoDeaths7DaysMAVG=0, NewConfirmedCasesPerMLagSlopeLM=quantile(obsSurvData$NewConfirmedCasesPerMLagSlopeLM, probs=0.25), TestPositivityRateLagSlopeLM=quantile(obsSurvData$TestPositivityRateLagSlopeLM, probs=0.25)",
         "DemGovernor = 0",
         "TrumpVote2016=quantile((obsSurvData$TrumpVote2016), 0.75)",
         "PctBlack=quantile((obsSurvData$PctBlack), 0.75)",
         "NewConfirmedCasesPerMLagSlopeLM=quantile(obsSurvData$NewConfirmedCasesPerMLagSlopeLM, probs=0.25)",
         "logNewDeathsPerM7DaysMAVG=log(quantile(obsSurvData$origNewDeathsPerM7DaysMAVG, probs=0.25)), NoDeaths7DaysMAVG=0",
         "TestPositivityRateLagSlopeLM=quantile(obsSurvData$TestPositivityRateLagSlopeLM, probs=0.25)",
         "PopDensity=quantile(obsSurvData$PopDensity, probs=0.25)"
         )


## Scenario names for coxed
nameScen <- c("Rep Gov + Trump Voters",
              "Falling Deaths, Pos & Cases",
              "Republican Governor",
              "More Trump Voters",
              "Higher Black Population",
              "Downward Trend New Cases",
              "Declining Rate of New Deaths",
              "Downward Trend Positivity",
              "Lower Population Density"
    )


## Make Figure 3 using coxed, which should then be polished in Adobe Illustrator.
## Note: the object ame contains the point estimates added to the polished plots
ame <- coxphAMEplot("Figure3", res1a, obsSurvData, obsSurvData$StateName, nameScen, cf2, cf1, ciMethod="studentized",
                    sims=1000,
log=FALSE,
                    pch=c(15, 16, 16, 16, 15, 17, 15, 16, 16),                    
                    col=c(nicered, niceblue, "black", niceblue, nicered, nicepurple, nicered, niceblue, niceblue),
                    limits=c(-0.35, 17.85),
                    at=c(0, 3.5, 7, 10.5, 14, 17.5),
                    labs=list(#c(" "),
                              c("0 days", "earlier"),
                              c(" "),
                              c("7 days", "earlier"),
                              c(" "),
                              c("14 days", "earlier"),
                              c(" ")
                              ),
                    title = "First Easing of Social Distancing Policy Expected...",
                    fancy=fancy
                    )


###################################################################################
## The rest of the code makes Figure 6, which shows many robustness checks


## Re-estimate baseline model with alternative sources of public health variables;
## compute counterfactuals, and save for plotting

## Create dataset using New York Times data
allsurvDataNYT <- mData
allsurvDataNYT <- allsurvDataNYT[(allsurvDataNYT$DaysSinceFCT>=startDate), ]
allsurvDataNYT <- allsurvDataNYT[(allsurvDataNYT$DaysSinceFCT<=endDate), ]
allsurvDataNYT <- allsurvDataNYT[(allsurvDataNYT$PolicyType=="GathRestrictRecomAny")|
                                 (allsurvDataNYT$PolicyType=="AnyBusinessClose")|
                                 (allsurvDataNYT$PolicyType=="RestaurantRestrict")|
                                 (allsurvDataNYT$PolicyType=="BarRestrict")|
                                 (allsurvDataNYT$PolicyType=="StayAtHome")
                                ,]
allsurvDataNYT <- allsurvDataNYT[allsurvDataNYT$COVIDSource=="NYT",]

## Create dataset using Covid Tracking Project data
allsurvDataCTP <- mData
allsurvDataCTP <- allsurvDataCTP[(allsurvDataCTP$DaysSinceFCT>=startDate), ]
allsurvDataCTP <- allsurvDataCTP[(allsurvDataCTP$DaysSinceFCT<=endDate), ]
allsurvDataCTP <- allsurvDataCTP[(allsurvDataCTP$PolicyType=="GathRestrictRecomAny")|
                                 (allsurvDataCTP$PolicyType=="AnyBusinessClose")|
                                 (allsurvDataCTP$PolicyType=="RestaurantRestrict")|
                                 (allsurvDataCTP$PolicyType=="BarRestrict")|
                                 (allsurvDataCTP$PolicyType=="StayAtHome")
                                ,]
allsurvDataCTP <- allsurvDataCTP[allsurvDataCTP$COVIDSource=="CTP",]

## Alternative datasets (replacing public health data source)
AltDataPH <- c("allsurvDataCTP", "allsurvDataNYT")

## Model formula for alternative public health data robustness checks
AltModelsPH <- c(sdEasingPH ~ DemGovernor + TrumpVote2016 + PctBlack + TestPositivityRateLagSlopeLM + logNewDeathsPerM7DaysMAVG +
                     NoDeaths7DaysMAVG + NewConfirmedCasesPerMLagSlopeLM + log(PopDensity) +
                     strata(PolicyType, SubstatePolicyDateSubNonReligIndoorEasing) + cluster(as.factor(StateName))
                 )

## Loop over possible alternative models (j) and datasets (i)
idx <- 0
resRobPH <- vector("list", length(AltDataPH)*length(AltModelsPH)) 
xprePH <- xpostPH <- rep(NA, length(AltDataPH)*length(AltModelsPH))
for (j in 1:length(AltModelsPH)) {
    for (i in 1:length(AltDataPH)) {

        sdEasingPH <-  Surv(time=get(AltDataPH[i])$DaysSinceFCT-1,
                           time2=get(AltDataPH[i])$DaysSinceFCT,
                           event=get(AltDataPH[i])$PolicySubNonReligIndoorEaseDV)

        idx <- idx+1
        resRobPH[[idx]] <-  coxph(AltModelsPH[[j]], data=get(AltDataPH[i]), x=TRUE, y=TRUE)
    

        ## "Before" scenario (factors theoretically disencouraging easing)
        xprePH[idx] <-  with(get(AltDataPH[i]),
                             c(log(quantile(origNewDeathsPerM7DaysMAVG, probs=0.75)),
                               log(quantile(origNewDeathsPerM7DaysMAVG, probs=0.75))
                               )[j]
                )

        ## "After" scenario (factors theoretically encouraging easing)
        xpostPH[idx] <- with(get(AltDataPH[i]),
                             c(log(quantile(origNewDeathsPerM7DaysMAVG, probs=0.25)),
                               log(quantile(origNewDeathsPerM7DaysMAVG, probs=0.25))
                               )[j]
                             )
    }
}


## Names of public health variable scenarios
robnamePH <- c(
    #"use JHU data for cases",
    "use CTP data for cases",
    "use NYT data for cases"
    )



## Perform remaining robustness checks


## Alternative samples:  leave out each policy area; leave out states that used substate policy
## First, create these datasets

## Leave out gatherings
allsurvData1 <- allsurvData[(allsurvData$PolicyType=="AnyBusinessClose")|
                            (allsurvData$PolicyType=="RestaurantRestrict")|
                            (allsurvData$PolicyType=="BarRestrict")|
                            (allsurvData$PolicyType=="StayAtHome")
   ,]
obsSurvData1 <- allsurvData1[!is.na(allsurvData1$PolicySubNonReligIndoorEaseDV),]
sdEasing1 <-  Surv(time=obsSurvData1$DaysSinceFCT-1,
                   time2=obsSurvData1$DaysSinceFCT,
                   event=obsSurvData1$PolicySubNonReligIndoorEaseDV)

## Leave out business closure
allsurvData2 <- allsurvData[(allsurvData$PolicyType=="GathRestrictRecomAny")|
                            (allsurvData$PolicyType=="RestaurantRestrict")|
                            (allsurvData$PolicyType=="BarRestrict")|
                            (allsurvData$PolicyType=="StayAtHome")
                           ,]
obsSurvData2 <- allsurvData2[!is.na(allsurvData2$PolicySubNonReligIndoorEaseDV),]
sdEasing2 <-  Surv(time=obsSurvData2$DaysSinceFCT-1,
                   time2=obsSurvData2$DaysSinceFCT,
                   event=obsSurvData2$PolicySubNonReligIndoorEaseDV)

## Leave out restaurant restrict
allsurvData3 <- allsurvData[(allsurvData$PolicyType=="GathRestrictRecomAny")|
                            (allsurvData$PolicyType=="AnyBusinessClose")|
                            (allsurvData$PolicyType=="BarRestrict")|
                            (allsurvData$PolicyType=="StayAtHome")
                           ,]
obsSurvData3 <- allsurvData3[!is.na(allsurvData3$PolicySubNonReligIndoorEaseDV),]
sdEasing3 <-  Surv(time=obsSurvData3$DaysSinceFCT-1,
                   time2=obsSurvData3$DaysSinceFCT,
                   event=obsSurvData3$PolicySubNonReligIndoorEaseDV)

## Leave out bar restrict
allsurvData4 <- allsurvData[(allsurvData$PolicyType=="GathRestrictRecomAny")|
                            (allsurvData$PolicyType=="AnyBusinessClose")|
                            (allsurvData$PolicyType=="RestaurantRestrict")|
                            (allsurvData$PolicyType=="StayAtHome")
                           ,]
obsSurvData4 <- allsurvData4[!is.na(allsurvData4$PolicySubNonReligIndoorEaseDV),]
sdEasing4 <-  Surv(time=obsSurvData4$DaysSinceFCT-1,
                   time2=obsSurvData4$DaysSinceFCT,
                   event=obsSurvData4$PolicySubNonReligIndoorEaseDV)

## Leave out stay-at-home
allsurvData5 <- allsurvData[(allsurvData$PolicyType=="GathRestrictRecomAny")|
                            (allsurvData$PolicyType=="AnyBusinessClose")|
                            (allsurvData$PolicyType=="RestaurantRestrict")|
                            (allsurvData$PolicyType=="BarRestrict")
                           ,]
obsSurvData5 <- allsurvData5[!is.na(allsurvData5$PolicySubNonReligIndoorEaseDV),]
sdEasing5 <-  Surv(time=obsSurvData5$DaysSinceFCT-1,
                   time2=obsSurvData5$DaysSinceFCT,
                   event=obsSurvData5$PolicySubNonReligIndoorEaseDV)

## Include as initial easing both religious easing and easing to takeaway/outdoor
allsurvData6 <- allsurvData[(allsurvData$PolicyType=="GathRestrictRecomAny")|
                            (allsurvData$PolicyType=="AnyBusinessClose")|
                            (allsurvData$PolicyType=="RestaurantRestrict")|
                            (allsurvData$PolicyType=="BarRestrict")|
                            (allsurvData$PolicyType=="StayAtHome")
                           ,]
obsSurvData6 <- allsurvData6[!is.na(allsurvData6$PolicyFirstSubEaseDV),]
sdEasing6 <-  Surv(time=obsSurvData6$DaysSinceFCT-1,
                   time2=obsSurvData6$DaysSinceFCT,
                   event=obsSurvData6$PolicyFirstSubEaseDV)

## Use alternative data focusing on easing in counties with 
## %Black>=state mean and total county pop at least state mean
allsurvData7 <- allsurvData[(allsurvData$PolicyType=="GathRestrictRecomAny")|
                            (allsurvData$PolicyType=="AnyBusinessClose")|
                            (allsurvData$PolicyType=="RestaurantRestrict")|
                            (allsurvData$PolicyType=="BarRestrict")|
                            (allsurvData$PolicyType=="StayAtHome")
                           ,]
obsSurvData7 <- allsurvData7[!is.na(allsurvData7$PolicySubNonReligIndoorEasingAtLeastOneBlackCtyDV),]
sdEasing7 <-  Surv(time=obsSurvData7$DaysSinceFCT-1,
                   time2=obsSurvData7$DaysSinceFCT,
                   event=obsSurvData7$PolicySubNonReligIndoorEasingAtLeastOneBlackCtyDV)


## Create model formulas for above 7 alternative datasets
AltModels <- c(sdEasing1 ~ DemGovernor + TrumpVote2016 + PctBlack + TestPositivityRateLagSlopeLM + logNewDeathsPerM7DaysMAVG +
                   NoDeaths7DaysMAVG + NewConfirmedCasesPerMLagSlopeLM + log(PopDensity) +
                   strata(PolicyType, SubstatePolicyDateSubNonReligIndoorEasing) + cluster(as.factor(StateName)),

               sdEasing2 ~ DemGovernor + TrumpVote2016 + PctBlack + TestPositivityRateLagSlopeLM + logNewDeathsPerM7DaysMAVG +
                   NoDeaths7DaysMAVG + NewConfirmedCasesPerMLagSlopeLM + log(PopDensity) +
                   strata(PolicyType, SubstatePolicyDateSubNonReligIndoorEasing) + cluster(as.factor(StateName)),
               
               sdEasing3 ~ DemGovernor + TrumpVote2016 + PctBlack + TestPositivityRateLagSlopeLM + logNewDeathsPerM7DaysMAVG +
                   NoDeaths7DaysMAVG + NewConfirmedCasesPerMLagSlopeLM + log(PopDensity) +
                   strata(PolicyType, SubstatePolicyDateSubNonReligIndoorEasing) + cluster(as.factor(StateName)),

               sdEasing4 ~ DemGovernor + TrumpVote2016 + PctBlack + TestPositivityRateLagSlopeLM + logNewDeathsPerM7DaysMAVG +
                   NoDeaths7DaysMAVG + NewConfirmedCasesPerMLagSlopeLM + log(PopDensity) +
                   strata(PolicyType, SubstatePolicyDateSubNonReligIndoorEasing) + cluster(as.factor(StateName)),

               sdEasing5 ~ DemGovernor + TrumpVote2016 + PctBlack + TestPositivityRateLagSlopeLM + logNewDeathsPerM7DaysMAVG +
                   NoDeaths7DaysMAVG + NewConfirmedCasesPerMLagSlopeLM + log(PopDensity) +
                   strata(PolicyType, SubstatePolicyDateSubNonReligIndoorEasing) + cluster(as.factor(StateName)),
               
               sdEasing6 ~ DemGovernor + TrumpVote2016 + PctBlack + TestPositivityRateLagSlopeLM + logNewDeathsPerM7DaysMAVG +
                   NoDeaths7DaysMAVG + NewConfirmedCasesPerMLagSlopeLM + log(PopDensity) +
                   strata(PolicyType, SubstatePolicyDateSub1stEasing) + cluster(as.factor(StateName)),

               sdEasing7 ~ DemGovernor + TrumpVote2016 + PctBlack + TestPositivityRateLagSlopeLM + logNewDeathsPerM7DaysMAVG +
                   NoDeaths7DaysMAVG + NewConfirmedCasesPerMLagSlopeLM + log(PopDensity) +
                   strata(PolicyType,  SubstatePolicyDateSubNonReligIndoorEasingAtLeastOneBlackCty) + cluster(as.factor(StateName))
               )

## Now create objects for sensitivity checks that add additional covariates
AddedVars <- c(
    "IdeoFordCiti",
    "ElectYear",
    "PctHispanic",
    "AtLeastCollege",
    "UrbanPct",
    "PctElderly",
    "InsuredUnempRate4WeeksMAVG",
    "log(GSPPerCapita)",
    "TourismEmpPct",
    "SalesTaxRevenueShare",
    "PovertyPct",
    "PctNeighSubNonReligIndoorEasing", 
    "PctAltSubNonReligIndoorEasing",
    "log(MeanNeighNewDeathsPerM7DaysMAVG)"
               )


## Create counterfactual scenarios for these checks

## "Before" scenario (factors theoretically discouraging easing)
xpre1robust <- with(obsSurvData,
                    c(NA,
                      rep(NA, length(xprePH)),
                      quantile(IdeoFordCiti, probs=0.75),       ## Ideology
                      0,
                      quantile(PctHispanic, probs=0.25), ## Hispanic
                      quantile(AtLeastCollege, probs=0.75),      ## AtLeastCollege
                      quantile(UrbanPct, probs=0.75),          ## Urban
                      quantile(PctElderly, probs=0.75),          ## PctElderly
                      quantile(InsuredUnempRate4WeeksMAVG, probs=0.25),
                      log(quantile(GSPPerCapita, probs=0.75)),
                      quantile(TourismEmpPct, probs=0.25),          ## Tourism
                      quantile(SalesTaxRevenueShare, probs=0.25),
                      quantile(PovertyPct, probs=0.25),
                      quantile(PctNeighSubNonReligIndoorEasing, probs=0.25), ##PctNeighSub1stEasing
                      quantile(PctAltSubNonReligIndoorEasing, probs=0.25), ## PctAltSub1stEasing
                      log(quantile(MeanNeighNewDeathsPerM7DaysMAVG, probs=0.75)),
                      NA,
                      NA,
                      NA,
                      NA,
                      NA,
                      NA,
                      NA,
                      NA
                      ))

    
## "After" scenario (factors theoretically encouraging easing)
xpost1robust <- with(obsSurvData,
                     c(NA,
                       rep(NA, length(xpostPH)),
                       quantile(IdeoFordCiti, probs=0.25),       ## Ideology
                       1,
                      quantile(PctHispanic, probs=0.75), ## Hispanic
                      quantile(AtLeastCollege, probs=0.25),      ## AtLeastCollege
                      quantile(UrbanPct, probs=0.25),          ## Urban
                      quantile(PctElderly, probs=0.25),          ## PctElderly
                      quantile(InsuredUnempRate4WeeksMAVG, probs=0.75),
                      log(quantile(GSPPerCapita, probs=0.25)),
                      quantile(TourismEmpPct, probs=0.75),          ## Tourism
                      quantile(PovertyPct, probs=0.75),
                      quantile(SalesTaxRevenueShare, probs=0.75),
                      quantile(PctNeighSubNonReligIndoorEasing, probs=0.75), ## PctNeighSub1stEasing
                      quantile(PctAltSubNonReligIndoorEasing, probs=0.75), ## PctAltSub1stEasing
                      log(quantile(MeanNeighNewDeathsPerM7DaysMAVG, probs=0.25)),
                      NA,
                      NA,
                      NA,
                      NA,
                      NA,
                      NA,
                      NA,
                      NA
                      ))



## Set up storage for robustness runs
resRob <- vector("list", 8 + length(resRobPH) + length(AddedVars) )#+ length(AltModels))

## Load baseline in slot 1
resRob[[1]] <- res1a

## Public health alternatives
for (i in 1:length(resRobPH)) {
    resRob[[1+i]] <- resRobPH[[i]]
}

## Added variable scenarios: run
for (i in 1:length(AddedVars)) {
    modelRob <- as.formula(paste(as.character(lhs(model1a)),
                                 "~", as.character(rhs(model1a))[2],
                                 "+", as.character(rhs(model1a))[3],
                                 "+", AddedVars[i]))
    resRob[[1 + i + length(resRobPH)]] <- coxph(modelRob, data=obsSurvData, x=TRUE, y=TRUE)
}

## Alternative outcome measures: run
resRob[[1 + length(resRobPH) + length(AddedVars) + 1]] <- coxph(AltModels[[1]], data=obsSurvData1, x=TRUE, y=TRUE)
resRob[[1 + length(resRobPH) + length(AddedVars) + 2]] <- coxph(AltModels[[2]], data=obsSurvData2, x=TRUE, y=TRUE)
resRob[[1 + length(resRobPH) + length(AddedVars) + 3]] <- coxph(AltModels[[3]], data=obsSurvData3, x=TRUE, y=TRUE)
resRob[[1 + length(resRobPH) + length(AddedVars) + 4]] <- coxph(AltModels[[4]], data=obsSurvData4, x=TRUE, y=TRUE)
resRob[[1 + length(resRobPH) + length(AddedVars) + 5]] <- coxph(AltModels[[5]], data=obsSurvData5, x=TRUE, y=TRUE)
resRob[[1 + length(resRobPH) + length(AddedVars) + 6]] <- coxph(AltModels[[6]], data=obsSurvData6, x=TRUE, y=TRUE)
resRob[[1 + length(resRobPH) + length(AddedVars) + 7]] <- coxph(AltModels[[7]], data=obsSurvData7, x=TRUE, y=TRUE)


## Names of scenarios for plot
robname <- c("Baseline Model",
             robnamePH,
             "Add: Citizen Ideology",
             "Add: Election Year",
             "Add: % Hispanic",
             "Add: % with a College Degree",
             "Add: % Urban Population",
             "Add: % at least 70 Years of Age",
             "Add: Unemployment Claims",
             "Add: State GDP pc",
             "Add: % Employment in Tourism",
             "Add: Sales Tax Dependence",
             "Add: Poverty Rate",
             "Add: % Neigh States with Policy",
             "Add: % of Peer States with Policy",
             "Add: N 7-day MA of New Death/M",
             "Exclude: Gatherings",
             "Exclude: Business Closures",
             "Exclude: Restaurant Restrict",
             "Exclude: Bar Restrict",
             "Exclude: Stay-At-Home Orders",
             "Alt: Include relig & outdoor",
             "Alt: Easing sizable Black pop"
             )

             
robxNum <- c(NA, rep(NA, length(resRobPH)), rep(9, 14), rep(NA, 7))

## Produce Figure 6 (all robustness checks).  
##Figure requires additional polishing in Adobe Illustrator
hrRob2 <- plotCovidCPH("Figure6",  res1=resRob,
                       scenname=robname,
                       xNum=list(robxNum,
                                 1,
                                 5,
                                 3),
                       xpre1=list(xpre1robust,
                                  1,
                                  log(quantile(obsSurvData$origNewDeathsPerM7DaysMAVG, probs=0.75)),
                                  quantile(obsSurvData$PctBlack, probs=0.25)),
                       xpost1=list(xpost1robust,
                                   0,
                                   log(quantile(obsSurvData$origNewDeathsPerM7DaysMAVG, probs=0.25)),
                                   quantile(obsSurvData$PctBlack, probs=0.75)),
                      limits=c(0.1, 3.1),
                      at=c(0.5, 1, 1.5, 2, 2.5, 3.0),
                      labsX=c("0.5", "1.0","1.5", "2.0", "2.5", "3.0"),
                      labsT=list(c("half", "as", "likely"),
                                 c("just", "as", "likely"),
                                 c("50%", "more", "likely"),
                                 c("100%", "more", "likely"),
                                 c("150%", "more", "likely"),
                                 c("200%", "more", "likely")),
                      sortBy=0, logaxis=FALSE, rev=TRUE,
                      titleT="The first move to ease social distancing is...",
                      titleX="Hazard Ratio",
                      entryheight=0, topOffset=0.02, topIncr=0.022,
                      pch=list(rep(16, length(robname)),
                               rep(15, length(robname)),
                               rep(16, length(robname)),
                               rep(17, length(robname))
                               ),
                      col=list("black", nicered, niceblue, nicepurple),
                      robustness=TRUE,
                      plotwidth=1.15, filewidth=14,
                      fancy=fancy
                      )




